      program bayescat5i
      implicit real*8 (a-h,l,o-z)
      parameter(nplay=101,ndec=12,ngrid=51)
c
      character*24 fdate
      dimension bgrid(ngrid),pz(ngrid),w0(ngrid),w(ngrid,ngrid)
      dimension ti(nplay,ngrid,ngrid),tl(ngrid,ngrid)
      dimension px(ndec,2),pxm(ngrid,ngrid,ndec,2)
      dimension is(nplay,ndec),xc10(ngrid,ngrid),ic10(ngrid,ngrid)
      dimension q10(nplay,4)
      common /mats/is
      common /pmatrx/pxm
c
      open(unit=11,file='../ra2011.data')   !with card sleeves on boards
      open(unit=12,file='cat10.data')
c
      read(11,*)((is(i,n),n=1,ndec),i=1,nplay)
      read(12,*)((xc10(i,j),j=1,ngrid),i=1,ngrid)
c
c-----------------------------------------------------------------
      write(*,'('' ***** bayescat5i ******'',20x,a24,
     &//"    Risk & Ambiquity Experiment")')fdate()
      write(*,'(/"  All 2011 Sessions"//)')
c
      do 5 i=1,ngrid
         if ((i.eq.1).or.(i.eq.ngrid)) then
            w0(i) = 0.5d0
         else
            w0(i) = 1.d0
         endif
    5 continue
c
      dz=1.d0/float(ngrid-1)
      do 20 i0=1,ngrid
         i=ngrid+1-i0
         if (i.lt.ngrid) then
            pz(i) = 0.5d0 + 0.5d0*dz*(i-1)
         else
            pz(i) = 0.999d0
         endif
         gam = 2.d0*dlog(pz(i)/(1.d0-pz(i)))
         do 15 j=1,ngrid
            beta = dz*(j-1)
            bgrid(j) = beta
            call fpx(gam,beta,px)
            do 10 k=1,ndec
               pxm(i,j,k,1) = px(k,1)
               pxm(i,j,k,2) = px(k,2)
   10       continue
            tl(i,j) = 0.d0
            w(i,j) = w0(i)*w0(j)
   15    continue
   20 continue
c
      do 50 id=1,nplay
         cumi = 0.d0
         do 35 i0=1,ngrid
            i=ngrid+1-i0
            do 30 j=1,ngrid
               call fam(id,i,j,t)
               ti(id,i,j) = t
               cumi = cumi + t*w(i,j)
   30       continue
   35    continue
         do 45 i0=1,ngrid
            i=ngrid+1-i0
            do 40 j=1,ngrid
               ti(id,i,j) = ti(id,i,j)/cumi
               tl(i,j) = tl(i,j) + ti(id,i,j)
   40       continue
   45    continue
   50 continue
c
      do 80 id=1,nplay
         cumi = 0.d0
         do 65 i0=1,ngrid
            i=ngrid+1-i0
            do 60 j=1,ngrid
               pij = tl(i,j) - ti(id,i,j)
               cumi = cumi + w(i,j)*ti(id,i,j)*pij
   60       continue
   65    continue
         do 75 i0=1,ngrid
            i=ngrid+1-i0
            do 70 j=1,ngrid
               pij = tl(i,j) - ti(id,i,j)
               pij = 1.d2*w(i,j)*ti(id,i,j)*pij/cumi
               k10 = xc10(i0,j)+1
               q10(id,k10) = q10(id,k10) + pij
   70       continue
   75    continue
   80 continue
c
      write(*,'(/"10% Categories:")')
      write(*,'(/"ID     AA      L0      Other    EUM"/)')
      do 110 id=1,nplay
         cumi = 0.d0
         do 105 k=1,4
            cumi = cumi + q10(id,k)
  105 continue
         write(*,'(i2,2x,5f8.3)')id,(q10(id,k),k=1,4),cumi
  110    continue
c
      end
c***********************************************************************
      subroutine fam(id,i,j,t)
      implicit real*8 (a-h,o-z)
      parameter (nplay=101,ndec=12,ngrid=51)
      dimension pxm(ngrid,ngrid,ndec,2)
      dimension is(nplay,ndec)
      common /mats/is
      common /pmatrx/pxm
c
      t=1.d0
      do 10 k=1,ndec
         m = is(id,k)
         t = t*pxm(i,j,k,m)
   10 continue
c
      return
      end
c***********************************************************************
      subroutine fpx(gam,beta,px)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12)
      dimension px(ndec,2)
c
      y1 = 10.d0
      y2 = 12.d0
      y3 = 15.d0
      y4 = 5.d0
      y5 = 6.d0
      y6 = 8.d0
      pjb = (1.d0+beta)/3.d0
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y1
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(1,1) = za/(za+zb)
      px(1,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y1
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(2,1) = za/(za+zb)
      px(2,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y2
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(3,1) = za/(za+zb)
      px(3,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y2
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(4,1) = za/(za+zb)
      px(4,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y3
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(5,1) = za/(za+zb)
      px(5,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y3
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(6,1) = za/(za+zb)
      px(6,2) = zb/(za+zb)
c
      za = y1*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(7,1) = za/(za+zb)
      px(7,2) = zb/(za+zb)
c
      za = y2*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(8,1) = za/(za+zb)
      px(8,2) = zb/(za+zb)
c
      za = y3*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(9,1) = za/(za+zb)
      px(9,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y4*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(10,1) = za/(za+zb)
      px(10,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y5*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(11,1) = za/(za+zb)
      px(11,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y6*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(12,1) = za/(za+zb)
      px(12,2) = zb/(za+zb)

c      do 10 k=1,ndec
c         write(*,'(i2,2g12.4)')k,(px(k,j),j=1,2)
c   10 continue
c
      return
      end

